        subroutine geomopt_car(NATOM,Num_ele,NORB,coor,ATOMCHG,baselable)
        implicit none
        integer :: i,j,k,l,icy,ind
        integer :: NATOM,Num_ele,NORB,maxcyopt,icyfinal
        integer :: ATOMCHG(NATOM)
        integer :: kbasst(NATOM),kbased(NATOM),kcenter(NATOM)
        integer :: mopttol,mgrad,moptalgo
        real*8,parameter :: a = 0.5291772083 !bohr= ai*a
        real*8  :: Eelec,Enuc,hh,rmsgrad,rmsdisp
        real*8  :: S(NORB,NORB),Hcore(NORB,NORB)
        real*8  :: U(NORB,NORB,NORB,NORB)
        real*8  :: fmaxtol,frmstol,dismaxtol,disrmstol !!! four Limit of convergence
        real*8  :: coor(NATOM,3),coortemp(NATOM,3),grad_rhf(3*NATOM)
        real*8  :: step(3*NATOM),coorznow(NATOM,3),mcoorznow(NATOM,3)
        real*8,allocatable :: coorstep(:,:,:),gradstep(:,:),record(:)
        real*8  :: hess(NATOM*3,NATOM*3),xold(NATOM*3),xnew(NATOM*3),xnow(NATOM*3) !!for quasi-Newton method
        character*30 baselable
        character*3 yn1,yn2,yn3,yn4
        character*10 destol
        !!!set parameter
        mopttol=1  !!!收敛限loose
        moptalgo=1 !!!梯度下降法
        mgrad=2    !!!数值梯度
        maxcyopt=500
        !!!end set 
        write(21,'(25X,A40)')'***** Geometry Optimization*****'
        write(21,'(25X,A40)')'use Cartesian coordinates '
        yn1='NO';yn2='NO';yn3='NO';yn4='NO'
        if(mopttol==1)then              !!  loose
            destol='Loose'
            fmaxtol=0.0025d0            !!  MAX Force
            frmstol=0.0016667d0         !!  RMS Force
            dismaxtol=0.01d0            !!  MAX Displacement
            disrmstol=0.0066667d0       !!  RMS Displacement
        else if(mopttol==2)then         !!  tight
            destol='Tight'
            fmaxtol=1.5d-5              !!  MAX Force
            frmstol= 1d-5               !!  RMS Force
            dismaxtol=6d-5              !!  MAX Displacement
            disrmstol=4d-5              !!  RMS Displacement
         else if(mopttol==3)then        !!  very tight
            destol='Very Tight'
            fmaxtol=1.5d-6              !!  MAX Force
            frmstol=1d-6                !!  RMS Force
            dismaxtol=6d-6              !!  MAX Displacement
            disrmstol=4d-6              !!  RMS Displacement
        end if
        write(21,*)'Geometry Optimization'
        write(21,*)'Max cycle is',maxcyopt
        write(21,*)'Convergence threshold value is',' ',destol
        allocate(coorstep(NATOM,3,maxcyopt+1),gradstep(NATOM*3,maxcyopt+1),record(5*maxcyopt+5))
        gradstep=0.0d0
        coorstep(:,:,1)=coor/a
        hh=0.001d0
        write(21,'(2(A6,3X),3(A17,3X))')'NUB','NZ','X','Y','Z'
        do i=1,NATOM
           write(21,'(2(I6,3X),3(F17.12,5x))')i,ATOMCHG(i),coor(i,:)
        end do
        call CalEnucre(NATOM,coorstep,ATOMCHG,Enuc) 
        !write(*,*)Enuc
        call INTXC_ALLOCATE(NATOM)
        call BASIS_INT(kbasst,kbased,kcenter,NATOM,NORB,coor,ATOMCHG,S,Hcore,U,baselable)
        call INTXC_DEALLOCATE()
        call RHF(S,Hcore,NORB,NATOM,Num_ele,U,Eelec)
        !write(*,*)Eelec+Enuc
        record(1)=Eelec+Enuc
        !stop
        record(2:5)=0.0d0
        do icy=1,maxcyopt
                write(21,*)''
                write(21,*)'======================================='
                write(21,*)'Geometry Optimization step',icy 
                write(21,'(20X,A40)')'Current Cartesian &
                                      coordinates(unit:A)'
                if(icy==1)write(21,*) 'Initial Cartesian coordinates'
                write(21,'(2(A6,3X),3(A17,3X))')'NUB','NZ','X','Y','Z'
                do i=1,NATOM
                write(21,'(2(I6,3X),3(F17.10,5x))')i,ATOMCHG(i),coorstep(i,:,icy)*a
                end do
                call CalEnucre(NATOM,coorstep(:,:,icy),ATOMCHG,Enuc) !!BOHR
                call INTXC_ALLOCATE(NATOM)
                call BASIS_INT(kbasst,kbased,kcenter,NATOM,NORB,coorstep(:,:,icy)*a,ATOMCHG,S,Hcore,U,baselable)
                call INTXC_DEALLOCATE()
                call RHF(S,Hcore,NORB,NATOM,Num_ele,U,Eelec)
               if(icy==1)then
                        write(21,*)'Initial E_ele:',Eelec
                        write(21,*)'Initial E_nuc:',Enuc
                else
                        write(21,*)' E_ele:',Eelec
                        write(21,*)' E_nuc:',Enuc
                endif
                        write(21,*)'Hartree Fock Energy   :',Eelec+Enuc
!!! calculata Gradient
                !write(*,*)'stop'
                if(mgrad==1)then     !!!解析梯度
                  call calgrad_rhf
                else if(mgrad==2)then!!!数值梯度
                  call Numgrad(NATOM,Num_ele,NORB,coorstep(:,:,icy),ATOMCHG,baselable,grad_rhf)
                else 
                  stop "Bad Gradient method"
                endif
                gradstep(:,icy)=grad_rhf
                coorstep(:,:,icy+1)=coorstep(:,:,icy)!!!coorstep(:,:,1)=coor
!!! Optimization algorithm
                if(moptalgo==1)then!!!梯度下降法 
                write(21,*)'Method: Gradient descent method'      
                write(21,*)'step :',hh
                step=-grad_rhf*hh
                do i=1,NATOM
                        do ind=1,3
                        coorstep(i,ind,icy+1)=coorstep(i,ind,icy)+step((i-1)*3+ind)
                        enddo
                enddo
                else if(moptalgo==2)then
                write(21,*)'Method: BFGS'
                
                endif
                !! update coordinates       
            coor=coorstep(:,:,icy+1)
            rmsgrad=sum(sqrt(grad_rhf**2.0d0))/sqrt(3.0d0*natom)
            rmsdisp=sum(sqrt(step**2.0d0))/sqrt(3.0d0*natom)
!!!             
            record((icy-1)*5+1)=Eelec+Enuc
            record((icy-1)*5+2)=maxval(-grad_rhf)
            record((icy-1)*5+3)=rmsgrad
            record((icy-1)*5+4)=abs(maxval(step))
            record(icy*5)=rmsdisp
            write(21,*)'deltaE',record((icy-1)*5+1)-record((icy-2)*5+1)
            write(21,'(2x,A12,A12,3x,A13)')'Item','Value','Convergen?'
            if(record((icy-1)*5+2)<fmaxtol )yn1='YES'
            write(21,'(2x,A12,F16.10,5x,A3)') 'maxgrad:',record((icy-1)*5+2),yn1
            if(record((icy-1)*5+3)<frmstol )yn2='YES'
            write(21,'(2x,A12,F16.10,5x,A3)') 'rms of grad:',rmsgrad,yn2
            if(record((icy-1)*5+4)<dismaxtol)yn3='YES'
            write(21,'(2x,A12,F16.10,5x,A3)') 'max disp:',abs(maxval(step)),yn3
            if(record(icy*5)<disrmstol)yn4='YES'
            write(21,'(2x,A12,F16.10,5x,A3)') 'rms of disp:',rmsdisp,yn4
            if(record((icy-1)*5+2)<fmaxtol .and.record((icy-1)*5+3)<frmstol &
              .and.record((icy-1)*5+4)<dismaxtol.and.record(icy*5)<disrmstol)then
                write(21,*)'After',icy,'steps,Optimization converged.'
                icyfinal=icy
                exit
            end if
            if(icy==maxcyopt)write(21,*) 'WARNING: NOT convergence!! '
        end do !!! icy


        write(21,*)' '
        write(21,'(30X,A40)')'*** Geometry Optimization Summary ***'
        write(21,*)' The convergence threshold of Maximum Force:',fmaxtol
        write(21,*)' The SCF convergence threshold of RMS Force:',frmstol
        write(21,*)' The convergence threshold of Maximum Displacement :',dismaxtol
        write(21,*)' The SCF convergence threshold of RMS Displacement :',disrmstol
        write(21,*) 'Final Cartesian coordinates'
        write(21,'(2(A6,3X),3(A17,3X))')'NUB','NZ','X','Y','Z'
        do i=1,NATOM
                write(21,'(2(I6,3X),3(F17.10,5x))')i,ATOMCHG(i),coorstep(i,:,icyfinal)*a
        end do
        write(21,*) 'Final Energy: ',record((icyfinal-1)*5+1),'a.u.' 
        deallocate(coorstep,gradstep,record)
        endsubroutine geomopt_car        


!!!======================================================================!
        subroutine Numgrad(NATOM,Num_ele,NORB,coor,ATOMCHG,baselable,grad_rhf)
        implicit none
        integer :: NORB,NATOM,Num_ele
        integer :: kbasst(NATOM),kbased(NATOM),kcenter(NATOM)
        integer :: i,ind,ATOMCHG(NATOM)
        real*8,parameter :: a = 0.5291772083 
        real*8  :: S(NORB,NORB),Hcore(NORB,NORB)
        real*8  :: U(NORB,NORB,NORB,NORB)
        real*8  :: h,E2,E1,Eelec,Enuc,grad_rhf(3*NATOM) 
        real*8  :: coor1(NATOM,3),coor2(NATOM,3),coor(NATOM,3)
        character*30 baselable
        !!! Numerical gradient
        write(21,*)'Numerical gradient'
        grad_rhf=0d0
        h=1.0d-5
        coor1=coor
        do i=1,natom
            do ind=1,3
            write(21,*)'atom number',i,'dimension: ',ind
            coor(i,ind)=coor1(i,ind)-h
            call CalEnucre(NATOM,coor,ATOMCHG,Enuc)
            CALL INTXC_ALLOCATE(NATOM)
            CALL BASIS_INT(kbasst,kbased,kcenter,NATOM,NORB,coor*a,ATOMCHG,S,Hcore,U,baselable)
            CALL INTXC_DEALLOCATE()    
            call RHF(S,Hcore,NORB,NATOM,Num_ele,U,Eelec)
            E1=Eelec+Enuc
            coor(i,ind)=coor1(i,ind)+h
            call CalEnucre(NATOM,coor,ATOMCHG,Enuc)
            !call CINT_DRIVER
            CALL INTXC_ALLOCATE(NATOM)
            CALL BASIS_INT(kbasst,kbased,kcenter,NATOM,NORB,coor*a,ATOMCHG,S,Hcore,U,baselable)
            CALL INTXC_DEALLOCATE()
            call RHF(S,Hcore,NORB,NATOM,Num_ele,U,Eelec)
            E2=Eelec+Enuc
            grad_rhf((i-1)*3+ind)=(E2-E1)/(2*h)
            end do !!! ind
        end do    !!! i 
        write(21,'(A60)')'Cartesian Gradient(Numerical gradient)&
                 dE/dX(Hartrees/Bohr)'
        write(21,'(2(A6,3X),3(A17,3X))')'NUB','NZ','X','Y','Z'
        do i=1,NATOM
           write(21,'(2(I6,3X),3(F17.12,5x))')i,ATOMCHG(i),grad_rhf((i-1)*3+1:(i-1)*3+3)
        end do  
        end subroutine   Numgrad
!!!======================================================================!
        subroutine Numhess()
        implicit none
        integer :: i,j,k
         
        endsubroutine


!!!======================================================================!
        subroutine calnucgrad(mgrad,NATOM,coor,grad_nuc)
        implicit none
        integer :: i,j,k
        integer :: mgrad
        integer :: NATOM,ind
        real*8  :: Raabb,tem
        real*8  :: coor(NATOM,3),ATOMCHG(NATOM)
        real*8  :: grad_nuc(3*NATOM)
!!! Gradient of nuclear repulsion energy
        if(mgrad==1)then
                write(21,*)'Gradient of nuclear repulsion energy'
        endif
        do i=1,NATOM
                do ind=1,3
                        tem=0.0d0
                        do j=1,NATOM
                        if(i==j) cycle
                        Raabb=(coor(i,1)-coor(j,1))**2+(coor(i,2)-coor(j,2))**2+&
                              (coor(i,3)-coor(j,3))**2
                        tem=tem + ATOMCHG(j)*((coor(j,ind)-coor(i,ind))/(Raabb)**1.5d0)
                        enddo
                enddo
                grad_nuc((i-1)*3+ind)=tem * ATOMCHG(i)
        enddo
        if(mgrad==1)then
        write(21,'(45X,A40)')'Cartesian Gradient dE/dX(Hartrees/Bohr)'
                write(21,'(2(A6,3X),3(A17,3X))')'NUB','NZ','X','Y','Z'
        do i=1,NATOM
           write(21,'(2(I6,3X),3(F17.12,5x))')i,ATOMCHG(i),grad_nuc((i-1)*3+1:(i-1)*3+3)
        end do
        endif
        endsubroutine calnucgrad
!!!======================================================================!
        subroutine CalEnucre(NATOM,coor,Zn,Enuc) !!!coor unit is BOHR
        implicit none
        integer :: i,j,NATOM
        integer :: Zn(NATOM)
        real*8  :: coor(NATOM,3)
        real*8  :: Raabb,Enuc
!!! Nuclear repulsion energy
        Enuc=0d0
        do i=1,NATOM
                do j=1,NATOM
                  if(i==j) cycle
                  Raabb=(coor(i,1)-coor(j,1))**2+(coor(i,2)-coor(j,2))**2+&
                        (coor(i,3)-coor(j,3))**2
                  Enuc=Enuc+0.5d0*Zn(i)*Zn(j)/sqrt(Raabb)
                 end do
        end do
        end subroutine   CalEnucre
 
